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Abstract Heterogeneity of neural attributes has recently gained a lot of attention and is increasing 
recognized as a crucial feature in neural processing. Despite its importance, this physiological feature 
has traditionally been neglected in theoretical studies of cortical neural networks. Thus, there is still a 
lot unknown about the consequences of cellular and circuit heterogeneity in spiking neural networks. 
In particular, combining network or synaptic heterogeneity and intrinsic heterogeneity has yet to be 
considered systematically despite the fact that both are known to exist and likely have significant roles 
in neural network dynamics. In a canonical recurrent spiking neural network model, we study how these 
two forms of heterogeneity lead to different distributions of excitatory firing rates. To analytically char¬ 
acterize how these types of heterogeneities affect the network, we employ a dimension reduction method 
that relies on a combination of Monte Carlo simulations and probability density function equations. 
We find that the relationship between intrinsic and network heterogeneity has a strong effect on the 
overall level of heterogeneity of the firing rates. Specifically, this relationship can lead to amplification 
or attenuation of firing rate heterogeneity, and these effects depend on whether the recurrent network is 
firing asynchronously or rhythmically firing. These observations are captured with the aforementioned 
reduction method, and furthermore simpler analytic descriptions based on this dimension reduction 
method are developed. The final analytic descriptions provide compact and descriptive formulas for how 
the relationship between intrinsic and network heterogeneity determines the firing rate heterogeneity 
dynamics in various settings. 

Keywords Leaky integrate-and-fire • Recurrent E/I Network • Intrinsic Heterogeneity • Network 
Heterogeneity • Dimension Reduction 


1 Introduction 

Theoretical studies of spiking neuronal networks have been extremely valuable for experimentalists and 
theoreticians. Uncovering the underlying mechanisms of complex phenomena in neural circuits often 
requires theory and/or computation. In this vein, this paper focusses on the effects of heterogeneous 
neural attributes in model neural networks. Heterogeneity is an undeniable physiological feature that has 
often been ignored in theoretical studies because it complicates theoretical analyses. Not only is there 
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ample experimental evidence for heterogeneity of neural networks at many scales (Markram et al 1997 


Parker 2003 Marder and Goaillard 2006 Bremaud et al 2007), but the importance of heterogeneity in 


neural computational is becoming more apparent. Indeed, a combination of theoretical and experimental 


studies on neural networks have demonstrated the value of heterogeneity (Ostojic et al, 2009 Hermann 


and Touboul 2012). In particular, theoretical studies have shown that heterogeneity can generally lead 


to efficient neural coding ( 

Shamir and Sompolinsky 

2006 Chelaru and Dragoi 

2008 Padmanabhan 

and Urban 2010 Marsat and Maler 2010 Mejias and Longtin 

2012 Tripathy et al 2013). However, 


there is still not a lot known about heterogeneity. Specifically, how do different sources of heterogeneity 
interact and lead to different neural network activity? 


We make a distinction between different sources of heterogeneity, addressing two forms: intrinsic 
and network heterogeneity, both of which are known to exist. Intrinsic heterogeneity are differences due 


to cellular properties that exist without coupling to other neurons (Marder and Goaillard 2006 Pad- 


manabhan and Urban 2010), for example the membrane time constant, threshold for spiking, reversal 


potentials, etc. Network heterogeneity, that is heterogeneity induced by coupling in a neural network, can 


arise from differences in synaptic coupling between neurons (Parker 2003 Marder and Goaillard 2006 


Bremaud et al 2007 Oswald et al 2009). To the best of our knowledge, the physiological relationship 


between these two sources of heterogeneity are not known. Therefore, we systematically study the effects 
of these forms of heterogeneity on a canonical recurrent spiking neural network. The main motivation 
for studying the relationship between different heterogeneous components is to provide a framework for 
possibly reconciling experimental measurements of multiple neural attributes; recent theoretical studies 
have shown that the network components can interact nonlinearly with surprising results ( [Marder and 


Goaillard 2006 Mejias and Longtin 2014 Hunsberger et al 20I4| ) (see Discussion section). The results 


of our study clearly show how multiple components effect the firing rate variability and might apply to 
experiments that measure the heterogeneity of these (or possibly other) neural attributes. 


The network we consider is noisy with variable spiking similar to that of real cortical neurons, and 
is excitable (i.e., neurons only fire with noise and/or synaptic coupling). We analyze how intrinsic and 
network heterogeneity together alter the dynamics of strongly coupled networks of neurons in various 
regimes ranging from asynchronous to rhythmic (i.e, ’ping’ network) firing. We focus on the dynamics of 
the excitatory neurons because they are the predominate cells for propagating signals to different layers 
in the cortex. Unsurprisingly, we find that when both intrinsic and network heterogeneity increase inde¬ 
pendently (i.e., when there is no relationship between them), the excitatory firing rates tend to have a 
larger range. However, for a fixed level of heterogeneity, the relationship or correlation between intrinsic 
and network heterogeneity strongly affects the overall range of excitatory firing rates. Moreover, these 
effects depend on what regime the neural network is in: during rhythmic firing, excitatory firing rate 
ranges decrease when intrinsic and network heterogeneity correlation increases; during asynchronous fir¬ 
ing, excitatory firing rate ranges increase when intrinsic and network heterogeneity correlation increases. 


To better understand these observations, we implement a dimension reduction method that relies on 
a combination of Monte Carlo simulations and analytic reductions. The reduction theory is based in part 


on our previous work (Ly 2013 Nicola et al 2015), and on the work of others (Moreno-Bote and Parga 


2006 Nesse et al 2008), where particular state variables are assumed to be slow and thus decoupled 


from other variables. Fortunately, the dimension reduction method captures the qualitative trend of the 
range of excitatory firing rates as heterogeneity is changed. This further inspires a simpler yet more 
revealing analytic description for how intrinsic and network heterogeneity combine to yield different 
ranges of excitatory firing rates. This study gives a more complete understanding of how heterogeneities 
interact and result in modulation of the firing rate statistics, which may ultimately lead to a better 
understanding of neural coding in coupled neural networks. 
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2 Methods 


2.1 Recurrent spiking LIF network 


The recurrent spiking network of excitatory (E) and inhibitory (I) neurons are modeled as leaky 
integrate-and-fire (LIF) neurons. The intrinsic and network heterogeneity are modeled simply by two 
parameters that are allowed to vary among the neurons. Other modeling studies impose heterogeneity 


in the response property (Tripathy et al 2013), e.g., orientation tuning (Shamir and Sompolinsky 2006 


Chelaru and Dragoi 2008). The models here effectively have heterogeneous response properties, but our 


focus is on two different sources that lead to that property. We model the intrinsic heterogeneity by 


having different voltage thresholds for spiking ( 

Mejias and Longtin 2012 Yim et al 

2013 

), equivalent to 

how many have incorporated intrinsic heterogeneity ( 

Strogatz and Mirollo 1991 

Chow 1998 Burton 


et al 2012 Ly et al 2012). The network heterogeneity is modeled by scaling the synaptic input by 
a value, effectively making each neuron receive different levels of network input. There is evidence in 


slice recordings that the probability of connection depends on distance (Oswald et al 2009 Levy and 
Reyes 2012), although we are not taking into account spatial dynamics, this model is plausible assuming 


synaptic strengths do not all inversely scale with connection probability. Moreover, there is abundant 
evidence for differences in synaptic coupling between neurons ( [Parker 2003[ [Marder and Goaillard 2006 


Bremaud et al 2007). 


The equations for the excitatory neurons indexed by j G {1, 2,..., A/'e} are: 
dv ■ 

= -Vj - gie{t){Vj - Si) - gee{t){Vj - Se) + 
(refractory period) ^ + ^ref) = 0 

drij 


= -Vj + 


dt 

9ee{t) = Qj 


gei{t) = 


i'efpresyn E cells} 

'lei 


Pei^i 


E 


Gk'{t) 


Td- 


/c'efpresyn I cells) 

Gj + Aj 


dG^ 

dt 

= -Aj+Tra^5{t-ti) 


( 1 ) 


where the leak, inhibitory and excitatory reversal potentials are 0, f/, and £e^ respectively with £i < 
< £e (the voltage is scaled to be dimensionless so that a leak/resting value of -65 mV maps to 0 and 
a threshold voltage of -55 mV maps to 1 (see Table 0)- ^j{t) are uncorrelated white noise processes, 
Pxy is the proportion of neuron type y (randomly chosen) that provides presynaptic input to neuron 
type X {x,y G {e,i}). The second line in the equations describes the refractory period at spike time 
t*: when the neuron’s voltage crosses threshold Oj (intrinsic heterogeneity), the neuron goes into a 
refractory period for Tref where the voltage is undefinecQ after which we set the neuron’s voltage to 
0. In the last line, ti denotes the spike times of the excitatory neuron. There are two factors in the 
equation for the total synaptic conductances {pee and Pei)'- qj and ; the latter does not depend 

on the individual neuron and is the same across the entire (E) population. However, qj introduces 
network heterogeneity by scaling both excitatory and inhibitory synaptic inputs. This form of network 
heterogeneity is loosely motivated by recent results by Xne et ^ ( |2Q14 ), who found that pyramidal 
neurons receive relatively similar proportions of excitation and inhibition in layer 2/3 of mammalian 
visual cortex (i.e., some cells receive more E/I while some cells receive less E/I). The qj factors are a 


In refractory, the other variables are governed by their ODEs 
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Table 1 Parameters for all simulations 


Parameter 

'T’m 

'Tref 

Si 

Se 

Tn 

Pxy 

For E and I: 

20 ms 

2 ms 

-0.5 

6.5 

5 ms 

0.2 


Parameter 

Td 

Tr 

a 

^efi 

E cells 

1 ms 

5 ms 

1 

800 

I cells 

2 ms 

10 ms 

2 

200 


Table 2 Parameters for specific regimes 


Regime: 

lei 

lie 

lee 

la 

(Te 

(Tl 

Noisy Rhythm 

10 

8 

12.25 

5 

2.5 

3 

Asynchronous 

10 

8 

0.05 

5 

3.5 

4 

Sharp Rhythm 

10 

8 

11.5 

5 

2.55 

- 


straight forward way to capture the different levels of ’balanced’ input (see equation (§) . The term 
network heterogeneity is often used to mean that the structure of the network is heterogeneous; here, 
we use that term to mean that the network activity induces heterogeneity via network inputs. 
Similarly, for the inhibitory neurons indexed by k G {1, 2,..., A^^}, the equations are: 


= -Vk - gii{t){vk -£l)- 9ei{t){vk - £e) + crirjkit) 
Vk{t*) > 1 (refractory period) => Vj{t* + Tref) = 0 

= -Vk + 


9ie{t) 


9ii{t) 




Tie 


Pie^e 

Tii 


E 


Gk> (t) 


Pii^i 


/c'efpresyn I cells} 


/c'etpresyn I cells} 


= -br/c + Ak 

at 

dAif . V—> . 

r^-— = -Ak + Tra^5{t - ti) 


( 2 ) 


Notable differences compared to excitatory neurons are that the threshold values are all equal to 1, 
and there is not a factor that scales the presynaptic inputs from the network. Although one could 
in principle relax these assumptions and augment the subsequent theory in a standard way, we made 
this choice because the results in this paper would not be diminished, and to avoid distracting from our 
focus on excitatory neuron behavior. The parameter values for all of the figures are in Table 

We consider two regimes of this model: (i) noisy rhythm, where the power spectrum is larger 
for certain frequency values, and (ii) asynchronous that has a flat power spectrum. Figure shows a 
comparison of various quantities in the three regimes considered in this paper (see Table for parameter 
values of different regimes). 


Monte Carlo Simulations: Monte Carlo simulations were run for 100 s of simulation time for ten realiza¬ 
tions (Fig. has only one realization) with a time step At = 0.2 ms, and the firing rates statistics were 
binned in non-overlapping 1 ms time windows. We use a common estimate of the standard deviation of 
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the firing rate calculations with Monte Carlo simulations across the n = 10 realizations: 


Kj) 




1 n 


i {j ) ^rr 


m=l 


the gray shaded regions in panels A and C of Figures |2}|^ represent 1 standard deviation above and 
below the sample mean: ± ^u{j) • The Monte Carlo simulation plots in panel B of Figures 2144 did 

not include standard deviation regions because the plots would be harder to see. The error bars of the 
Monte Carlo simulations in panel D of Figures [2]-[^ and Figure [^, C, E, represent an estimate of the 
standard deviation of: max u — min u. We use the following estimate: 


(max) — zy (min) 


' ^/(max) 


^^/(min) 


even though cr^(max)-zy(min) = Ycr^(max)^ +cr^(min)^ “ 2 (^ 01 ;(z/(max), z/(min)); the reason for this is be¬ 
cause estimating 2Co'z;(z^(max), z^(min)) would require storing an additional 0{N^) sized vector, as well 
as O(n^) simulations for similar order accuracy (for every single parameter set), and the computation 
times are already quite long. 


2.2 Model with a sharper rhythm 


Another model considered is one with less variable inhibitory firing that ultimately leads to sharper 
rhythms in the excitatory neurons. The reason for an observed sharper rhythm (Fig.[^) is that in these 
recurrent networks, the inhibitory neuron firing silence and thus shape the rhythm of the excitatory 


neurons (Borgers and Kopell 2003 Economo and White, 2012). The only change is in the inhibitory 


neuron’s voltage equation which no longer has a noise term, but rather a deterministic drift to a sub¬ 
threshold target voltage Sdet- 

'^m ~ ^/) ^e) Odeti^k ^det) 

we set £det = 0.9 and gdet = 2. The regime considered has a strong and relatively regular oscillation 
(Eigure [fright column sharp rhythm). 

In addition to showing distinct characteristics compared to the other regimes (noisy rhythm, asyn¬ 
chronous), the main motivation for considering this additional model is that such recurrent networks 
with a sharp gamma rhythm are commonly studied and known to be important for coding in many areas 


2012 ). 


of the cortex (Borgers and Kopell 2003 Wang 2010 Buzsaki and Wang 2012 Economo and White 


The following two subsections describe the way both network and intrinsic heterogeneity are modu¬ 
lated in this paper. 


2.3 Changing the level of intrinsic and network heterogeneity independently 

The two heterogeneous parameters are varied to yield significant changes in the range of firing 

rates. The means of both q and 6 are set to 1, and the parameters G [0,1] and ae G [0,1] quantify 
the level of the network and intrinsic heterogeneities, respectively, in the following way: 

q ^ 1 + CTg * (Z// — 0.5) 

0 ^ 1 (j0 ^ J\f 
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( 3 ) 

( 4 ) 






































where U is the uniform distribution on [0,1], and A/’ is a truncatecQ normal distribution with mean 0 
and standard deviation 0.08. When chosen independently, the correlation between these two vectors will 
be small and theoretically zero. 


2.4 Changing the correlation between intrinsic and network heterogeneity 

We consider another way to change the heterogeneity where the overall level is approximately the same 
but the correlation between q and Q is set to a prescribed value. Given the vectors q and we fix q 
to the same values but transform d so that the Pearson’s correlation coefficient is q G ( — 1,1) in such 
a way that the transformed vector has the same mean and variance as d. The details for how this is 
accomplished are described in the Appendix. 


2.5 High dimensional probability density equation 


The recurrent coupled stochastic network in section |2.1| is difficult to describe theoretically. A common 
method uses probability density functions (p.d.f.), or a population density methods, where the proba¬ 
bility of a neuron being in a particular state has a corresponding equation. The variables in the popu¬ 
lations are no longer tracked individually, but rather captured by a p.d.f.; for example, {Vj^Gj, Aj^rjj) 
for j = l,2,...,A^e are captured with a function of {ve^ QEiCiE^VE)- The two forms of heterogeneity 
introduce even more dimensions than the usual state variables. For simplicity, one can track a family 
of probability density functions for each (qj^Oj) pair or each distinct neuron. The subsequent equations 
are a good approximation to the coupled network §-(§ with the following assumptions: 

(i) finite size effects are negligible (A'e/i ^ 1) 

(ii) the firing rate of presynaptic neurons is governed by a Poisson process 

(iii) the population firing rate averaged over q and 0 is a good approximation to the average presynaptic 
input rate 

(iv) a single p.d.f. function is adequate to describe the population behavior, and the heterogeneity is 
driven by {qj^Oj) 


The first two assumptions are standard in this framework, while the last two assumptions has been suc¬ 
cessfully used (Ly 2014), where a family of probability density functions were indexed by the quenched 
heterogeneity values. Even though these assumptions are violated, the following equations are key for 
the reduced descriptions in sections |2.6[ |3.2f[3.3| 

For each pair of values {qj^Oj), the probability density function p is defined by: 


J piVE,VfE,Vl,Vfi,t)dVEdWEdVidwi=PT(^[vE{t),VfE{t),Vi{t),Wi{t)) e J?) 

where wx denotes the other states variables of the corresponding neuron type X G {E, /}, consisting of 
conductance, colored noise: wx = {gx^ctx^Vx)- The evolution of the p.d.f.’s is governed by a continuity 
equation and boundary conditions: 


^ The middle 98.76% is included, so for ag = 1, 6 G [0.8,1.2] 
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dp 

dt 

J 


^ = -V-J 


•— {Jvet Jg El El JrjET Jgn Ja 11 Jrii) 


JvE •=- [^E + qieidlivE — £l) + qjeeSEivE — £e) + <^EVe] P 


Jdt • 


J. 


gx 


(5) 

(6) 

(7) 

(8) 
(9) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 

The definitions of qxy in the LIF equations Q-© result in a total conductance of ^xy9y on average. 
Note that with a refractory period Tref > 0, the system of equations should also include a refractory 


-b/ + ~ ^l) + — £e) + CFi 9 i] P 

'T~rn 

ax\p 

pax 

it) / 

J a\ 


--[gx 

Ed 

ax , , 

-h t^x{ 

Tr 


p{. ■., a’x, ■ ■ ■) da'x 


Jr)x ■= - Vxp- 


ax—ax 

1 d^p 


Exit) ;= — 


UII 


Tn dr]‘ 


Jvx = 6>, Wx) dwx dq dO 


Jvxi^X = O.Wx^t) 
Jwx l^wx 


= Jvxi^X = 0,Wx,t+ r^e/) 


probability density that we do not state here (see the work of Tranchina and colleagues Nykamp and 
[Tranchina (2001); Haskell et al (2001); Apfaltrer et al (2006); Ly and Tranchina (2009) for further 
details). 


2.6 Reduction theory to describe firing rate dynamics 


We describe an insightful analytic reduction that captures how the range of excitatory firing rates change 
in different regimes. We focus on only the excitatory neurons, which have fewer state variables if the 
inhibitory population is ignored or assumed to be known. The problem with using the full p.d.f equations 
(§-([14]) is that the state variables are coupled, so we will formally assume that all of the state variables 
are known (given) except ve^ and solve for the steady-state firing rate as a function of the other state 
variables ( |Ly[|2Q13{ Nicola et al 2015). Note that other authors have employed a similar approach using 


an adiabatic or slow variable approximation in the context of stochastic spiking neurons (Moreno-Bote 
and Parga 2QQ6[ [Nesse et al 2008); very recently Hertag et al ( 2014| ) used this approach formally (see 


their equation (25)). 

We denote r as the approximation to the excitatory firing rate(s) ee- Assuming the other 
state variables are simply parameters, the deterministic firing rate of the equation 


'^m 


dVE 

dt 


= -ve - qgi{vE - £i) - q9~E{vE - £e) + 9e 


is straightforward to calculate, and given by 


ro{q,0;wE) 


0 , 


t log(^ 


l+q(g~E+g~i) 


y 


_ (li9E^E+9l£l) + VE _ 

9i9~E^ E+9~I £ l) + 'n~E ~^i^ + 9i9~E+9~l)) 


•r q{g~E£E+gi£i)+r]~E 
lYq{gYYgi) 

•r q{g~E£E+gi£i)+rfE 
lYq{g~E+gi) 


<0 

>0 


(15) 


The argument of the left-hand side is written in this way because we assume the q and 0 values are 
the primary sources of heterogeneity, rather than the external noise, finite size effects, random 
connectivity, etc. For exposition, we have absorbed the parameters and defined new variables with tildes: 
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g~E = leeOE^ Qi = leidE Ve = cte^e- The approximation (15) ignores the refractory period, which is 
accounted for via a transformation rather than using the refractory probability density. The inverse of 
the firing rate is the time between spikes, so the refractory period can be added to the time between 
spikes to yield: ro/(l + T^e/^o) as a simple approximation to the firing rate. Finally, the given state 
variables are integrated against their marginal density to get: 


r{q,e)=^ 


ro 


1+roTrefj J 1 + roTref 


/ 


ro 


-K9e,9i,Ve)<1we 


(16) 


There is a slight abuse of notation because the auxiliary variables ax effect the conductances but are 
not written in the previous equation; the emphasis is on how {qe^ gi^^E) directly effects r. Since the 
external noise is applied indiscriminately, 77 k is independent of the other variables and the marginal 
density factors into: 

K9~e,9i,Ve) =p{9'e,9i) -^—• 

cteve 


However, p{g~E^gi) is still not analytically tractable, leading us to rely on Monte Carlo simulations to 
numerically estimate p{g~E^gi)‘ 

The reduction method © -(p^ was implemented by relying on Monte Carlo simulations for p{g^E^gi)i 
and using the same vectors (q, 0) in the LIF simulations. Since p{g~E^gi) is the same for a given param¬ 
eter set, the range of firing rates is theoretically captured by the different (g^, Oj) values in equation ( [T^ 
(see blue curves in Figs.|2||^. 


3 Results 


We consider a recurrently coupled stochastic spiking neural network. Such networks have been ubiquitous 
in contemporary theoretical investigations. The class of networks considered here are widely used even 
though we do not include plasticity or detailed biophysical properties with different time scales. We are 
interested in how the firing rate of the excitatory population changes as the level of heterogeneity is 
varied in different regimes. Excitatory neurons are the focus because they are the predominate cells for 
propagating signals to different layers in the cortex. 

Figure highlights the different behaviors in three regimes considered. Representative raster plots 
of spikes are shown in panels A-C for both excitatory and inhibitory cells, incorporating both forms of 
heterogeneity with (q, 0) chosen independently and with = 1. In panels D-F, the power spectrum 

of the excitatory population firing rate for both heterogeneous (black) and homogeneous parameters 
(magenta) are shown; the thinner lines are the power spectrums of the individual excitatory neurons. 
The inset shows the autocorrelation function of the (E) population firing rate. The autocorrelation 
function of a stochastic process R{t) (e.g., excitatory population firing rate) is: 

A{t) = [R{r)R{t - r)] - E^[R(r)]^ (17) 


and the power spectrum is: 


P(co’) = / A{t)e 


I 


— i2'KUjt 


dt 


(18) 


These quantities illustrate the different dynamics of each neural network. Notice that the power spectrum 
of the individual neurons is consistently larger in value than the power spectrum of the population 
firing rate. This is because the population firing rate is averaged (smoothed) so that P{uj) —> 0 for 
large frequencies, whereas the spike train of the individual neurons is not smooth, consisting of O’s 
or S's ~ l/{At) that will always yield Pj > 0 as long as u; < l/{At) (the numerical limit because 
of the discretization). Specifically, as cj ^ l/(Z\t), Pj measures the average power of the spike train 
in a single time bin and thus converges to the firing rate of the individual neuron. The bottom row 
(panels G-I) show the distribution of the excitatory firing rates for each individual neuron, averaged 










Table 3 Firing Rate Values in Figure |i|G-I. The inhibitory firing rates are not shown in Figure [fj 


Regime: 


Mean u (Hz) 

[^miriT^max] (H^l) 

Range of u (Hz) 

Noisy Rhythm 

(Heterog.): E cells 

12.8 

[6.6, 23] 

16.4 

Noisy Rhythm 

(Homog.): E cells 

11.6 

[10.3, 12.7] 

2.5 

Asynchronous 

(Heterog.): E cells 

7.7 

[2.7, 14.2] 

11.4 

Asynchronous 

(Homog.): E cells 

7.4 

[6.5, 8.2] 

1.7 

Sharp Rhythm 

(Heterog.): E cells 

34.1 

[14.2, 64] 

49.8 

Sharp Rhythm 

(Homog.): E cells 

33.6 

[32.1, 35.5] 

3.5 

Noisy Rhythm 

(Heterog.): I cells 

17 

[15.9, 18.7] 

2.8 

Noisy Rhythm 

(Homog.): I cells 

16.1 

[14.8, 17.4] 

2.6 

Asynchronous 

(Heterog.): I cells 

19.7 

[18.5, 21.2] 

2.7 

Asynchronous : 

Homog. I cells 

19.6 

[17.9, 20.9] 

3 

Sharp Rhythm 

(Heterog.): I cells 

26.8 

[17.4, 32.4] 

15 

Sharp Rhythm 

(Homog.): I cells 

26 

[24.3, 27.9] 

3.6 


over time (simulations in Figure were performed for 100 s). The heterogeneous population naturally 
has a wider distribution compared to the homogeneous regime. At the population level, there are only 
minor differences between the homogeneous and heterogeneous regimes; indeed, the average firing rates 
(not the overall distribution), power spectrums, and autocorrelation functions are very similar. Thus 
enabling a systematic assessment of how intrinsic and network heterogeneity effect the spiking network, 
avoiding the complication of regime changes due to heterogeneity. Although there has been interesting 
work showing how heterogeneity can induce rhythms from asynchrony (i.e., bifurcations) ( [Hermann and 


Touboul 2012 Mejias and Longtin 2012), we do not directly address such dynamics here. Our study 


focuses on comparing the firing rate heterogeneity modulation within specific regimes. 

The neurons in all regimes considered are all excitable and receive external colored noise, resulting 
in irregular spiking. A common measure of variability is the Fano factor, defined as the variance of spike 
counts divided by the mean of spike counts in a time window. For all of the networks considered in this 
paper (e.g.. Fig. [^, the Fano factor of the E population spike counts is often greater than 1 and is at 
least 0.8 for heterogeneous and homogeneous networks, across all regimes, and for time windows ranging 
from 2 to 50ms (not shown). 


3.1 PDF framework captures firing rate modulation with heterogeneity 

The level of intrinsic and network heterogeneity were modulated in the recurrently coupled networks 
in the two ways previously described: i) independently choose the vectors (q, 6 ) with a prescribed level 
of heterogeneity determined by and ae (see and ii) for fixed values of change the 

correlation between (q, 0). In the noisy rhythm regime. Figure]^ shows the minimum and maximum 
excitatory firing rates for fixed values of = 0 and 1, while ao varies between 0 and 1 (black curves; 
dashed and solid). Not surprisingly, as ao (or cr^) increased, so does the range of (excitatory) firing 
rates with the minimum decreasing and maximum increasing. Note that we chose to plot the curves 
for fixed cr^; the results also hold with fixed ae and aq on the x-axis (not shown). Figure]^ shows the 
range of the firing rates (maximum minus the minimum) rather than the raw firing rate values; it is 
apparent that more intrinsic and/or network heterogeneity leads to more firing rate heterogeneity. The 
reduction theory described in section [2^ (based on both probability density equations and Monte Carlo 
simulations), and in particular equations for the approximation of the firing rates, are shown 

in the blue colored curves. In Figure [^, the theory does not provide a good quantitative match. This 
could be due to a variety of reasons: the 4 assumptions in section [23] are violated, the reduction method 
is known to be inaccurate compared to both the full PDF and Monte Carlo simulations. Fortunately, the 
reduction theory is able to capture the increase in the firing rate range in Figure]^, where aq^e ^ [0, !]• 
This result is indeed fortuitous given the inaccuracies of the PDF theory in capturing the raw firing 
rates. 
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Fig. 1 The three regimes considered are: noisy rhythm (left column, A, D, G), asynchronous (middle column, B, E, H), 
sharp rhythm (right column, C, F, I). The top row (A-C) has representative raster plots of spikes for both excitatory 
(black dots) and inhibitory (green dots) neurons with both intrinsic and network heterogeneity, showing distinct behavior 
depending on the regime. The middle row (D-F) shows the power spectrum of the excitatory population firing rate 
with both forms of heterogeneity (black) and homogeneous parameters (magenta); the power spectrum of the individual 
excitatory neurons are shown with thinner lines. The inset of each panel shows the autocorrelation function of the excitatory 
population rate. The bottom row (G-I) has histograms of the average firing rate for each excitatory neuron, with the 
heterogeneous network naturally having a wider distribution. The mean firing rate, minimum and maximum firing rates, 
and the range of the firing rates are displayed in TableThe intrinsic and network heterogeneity parameters were selected 
independently with aq = ag = 1 (see §-@)- The simulations were performed for a single realization of 100 s. 


Figure]^ shows the minimum and maximum firing rate with ample heterogeneity {aq = = 1, 

the most we considered), but the correlation between q and 9 (p(q, 0)) varied between (-1,1)- The 
comparison of the reduction theory (15)- (blue curves) and the simulations (black curves) is not 
accurate (as in Fig. [^), but does qualitatively capture the trend in the range of the firing rates (Fig. 
[^). The range of firing rates tends to decrease as p(q, 0) increases. Note that the range of firing rates 
can be very large and very small depending on p, even though aq = ao = 1; in fact, the ranges of firing 
rates are comparable to varying the overall level of heterogeneity: aq/Q G [0,1] (Fig. j^-B). In other 
words, a particular range of excitatory firing rates can arise from different levels of intrinsic and network 
heterogeneity, depending on their relationship. 
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Correlation of (0,q) Correlation of (0,q) 

Fig. 2 Noisy Rhythm regime: excitatory firing rates modulation with changes in intrinsic and network heterogeneity. 
A)-B) changing the level of intrinsic gq and network Gq heterogeneity by independently drawing (q, 0) (see §-@). 
A) The minimum and maximum firing rates are shown for the homogeneous network (cTg = cr^ = 0, dash curves) and 
completely heterogeneous (cTg = gq = solid curves) network. The simulation curves (black curves, dash and solid) 
are Monte Carlo simulations of equations (shaded regions denote standard deviations ). H ere, the theory uses a 

combination of dimensionally reduced PDF functions and Monte Carlo simulations (see ( |15| >- |T^ ). Although the theory 
does not quantitatively match the extreme values of the firing rates (A), it captures the trend of the firing rate range 
(maximum minus minimum) in panel B); there, each curve is for a fixed value of Gq ranging from 0 to 1, as cr^ ranges 
between 0 and 1. Unsurprisingly, as heterogeneity increases so does the firing rate range. The firing rates and the range 
vary appreciably over an order of magnitude. C)-D) changing the correlation q between (q, 0) with Gq = gq = 1. C) 
the theory (blue curve) does not quantitatively capture the actual firing rates as q varies between (-1,1), but they are 
comparable. D) the theory (solid) captures the trend in the simulated range of the firing rates (dots) as q varies between 
(-1,1). As Q increases, the range of firing rates tends to decrease. Gray shaded regions in A and C are an estimate of 
the standard deviation about the sample mean of the Monte Carlo simulations (100 s simulation for each realization, 10 
realizations total); er ror bars in D are estimates of the standard deviation about the sample mean of the range (see last 
paragraph of Section [ tT] for details). Shaded regions are omitted in B for readability. 


Similar comparisons are made for the two other regimes in Figures and In the asynchronous 
regime, as the the heterogeneity parameters are selected independently (Fig. Ih -B), we again see that 
more heterogeneity leads to a wider range of firing rates. Figure]^ shows = 0 and 1 split into two 
panels so it is easier to compare the theory (15)-(16) and simulations 0 (§• The quantitative match 
is not good (Fig. [^) as expected given the previous figure, but the trend is captured (Fig. |^, where 
(Tq^e ^ [O 5 !])• For fFe bottom row, we fix (jg = = 1 and let the correlation g between intrinsic and 

network heterogeneity vary between (-1,1)- Notice that the range of firing rates changes in a different way 
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Fig. 3 Asynchronous regime: excitatory firing rates modulation with changes in intrinsic and network heterogeneity 
[similar to Figure]^ A)-B) changing the level of intrinsic ag and network aq heterogeneity by independently drawing (q, 6 ) 
(see ([^-(^). A) The minimum and maximum firing rates are shown for the completely heterogeneous (aq = ag = solid 
curved network in the top panel and the homogeneous network (aq = ag = 0, dash curves) in the bottom panel. The 
simulation curves (black curves, dash and solid) are Monte Carlo simulations of equations ([T] >-( [^. Here, the theory uses 
a combination of dimensionally reduced PDF functions and Monte Carlo simulations (see ( |l^ - |l^ ). Although the theory 
does not quantitatively match the extreme values of the firing rates (A), it captures the trend of the firing rate range 
(maximum minus minimum) in panel B); there, each curve is for a fixed value of aq ranging from 0 to 1, as cr^ ranges 
between 0 and 1. Unsurprisingly, as heterogeneity increases so does the firing rate range. C)-D) changing the correlation 
g between (q, 6 ) with aq = ag = 1. C) the theory (blue curve) does not quantitatively capture the actual firing rates as g 
varies between (-1,1). D) the theory (solid) captures the trend in the simulated range of the firing rates (dots) as g varies 
between (-1,1). As g increases, the range of firing rates tends to increase. Cray shaded regions in A and C are an estimate 
of the standard deviation about the sample mean of the Monte Carlo simulations (100 s simulation for each realization, 10 
realizations total); er ror bars in D are estimates of the standard deviation about the sample mean of the range (see last 
paragraph of Section [ tT] for details). Shaded regions are omitted in B for readability. 


(Fig.l^). Here, as g increases, the range of firing rates increases in contrast to before where it decreased 
(Fig.l^). Moreover, we see the range of firing rates change by a factor of ^3, which interestingly is 
comparable to the firing rate range values when varying (jq/g independently (Fig. 

In Figure the sharp rhythm regime shows similar characteristics to Figure except for the 
following. The excitatory firing rate range is more sensitive to varying the heterogeneity parameters, 
and we see that the range of firing rates takes on much larger values. This is apparent both when aq and 
ag vary independently (Fig.|^-B) and when the correlation between q and 6 changes with aq = ag = 1 
(Fig. [^-D). The firing rate range changes by almost an order of magnitude. The other interesting 
thing about this regime is that the reduction theory (blue) matches the simulations much better than 
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Fig. 4 Sharp Rhythm regime: excitatory firing rates modulation with changes in intrinsic and network heterogeneity 
[similar to Figures A)-B) changing the level of intrinsic cfq and network cfq heterogeneity by independently drawing 
(q, 6) (see ^ -0) . i^The minimum and maximum firing rates are shown for the completely heterogeneous (cTg = ag = 1, 
solid curves) network in the top panel and the homogeneous network (aq = ag = 0, dash curves) in the bottom panel. The 
simulation curves (black curves, dash and solid) are Monte Carlo simulations of equations ([T| >-( ^. Here, the theory uses 
a combination of dimensionally reduced PDF functions and Monte Carlo simulations (see 1^). Although the theory 

does not quantitatively match the extreme values of the firing rates (A), it captures the trend of the firing rate range 
(maximum minus minimum) in panel B); there, each curve is for a fixed value of aq ranging from 0 to 1, as cr^ ranges 
between 0 and 1. Unsurprisingly, as heterogeneity increases so does the firing rate range. C)-D) changing the correlation 
g between (q, 0) with aq = ag = 1. C) again, the theory (blue curve) does not quantitatively capture the actual firing 
rates as g varies between (-1,1). D) the theory (solid) captures the trend in the simulated range of the firing rates (dots) 
as g varies between (-1,1). As g increases, the range of firing rates tends to decrease. Gray shaded regions in A and C are 
an estimate of the standard deviation about the sample mean of the Monte Carlo simulations (100 s simulation for each 
realization, 10 realizations total); er ror b ars in D are estimates of the standard deviation about the sample mean of the 
range (see last paragraph of Section [2.1 [ for details). Shaded regions are omitted in B for readability. 


the other two regimes (Fig. The match is particularly good for the maximum firing rate. Similar 

to the noisy rhythm regime, we see that as g increases, the range of firing rates decreases dramatically 
(Fig. 1^). One striking observation is that the PDF theory captures the firing rates much better in the 
sharp rhythm regime (Fig. EP -D) than the other two regimes (Fig. [^-D, Fig. |^-D). Although the 
underlying reason for this is difficult to determine exactly, a plausible explanation is that the overall 
higher firing rates and less noise in the sharp rhythm regime are consistent with the assumptions of the 
approximation in equation (15). 

The observation that more intrinsic and (uncorrelated) network heterogeneity leads to a wider range 
of firing rates (Figs. |2||4j panels A and B) is expected and does not require further analytical insight. 
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However, changing the correlation between q and 6 for fixed values of aq and cfq results in enlarged or 
diminished ranges of firing rates, depending on the regime (Figs.[2]-[^ panels C and D). How can these 
observations be reconciled? The next two sections provide further analysis to deeply understand 
these phenomena. 


3.2 Analytic description of heterogeneous firing rate range in rhythmic networks 

For many regimes and the different types of heterogeneity, the reduction method in section |2.6| does 
qualitatively capture the modulation of the range of firing rates, as shown in Figures |2}|^ panels B and 
D. This observation is the motivation for further analysis of equations that ultimately yields 

simple analytic formulas to account for how the firing rate ranges are effected by the relationship between 
intrinsic and networks heterogeneity. In sections |3.2| and |3.3[ we use the following variable substitutions 
to facilitate exposition of the analysis: 


xo:=g~E^gi (19) 

xi := g~ESE + gi^i ( 20 ) 


When the coupling parameters yield rhythmic firing (i.e., power spectrum of the population firing 
rate is not fiat), the net synaptic input is large on average (averaged over time and across excitatory 
neurons) and is much larger than when the network is in an asynchronous regime. Thus, we consider 
the large firing rate limit in the reduced theoretical description ([T^-(p!^. Furthermore, we ignore the 
effects of the refractory period r^e/ external noise r]~E^ and focus on the formula for the deterministic 


firing rate (15): 


^ro{q,0) 


1 + qxo 


log ( 


qxi-6{l-\-qxo) I 


qxi 


( 21 ) 


We assume the random state variables are parameters just like in the reduced description, which will 
enable us to focus on (g, 0) and determine how these two parameters effect the firing rate range. The 
two vectors (q, 0) are the main source of the firing rate heterogeneity. The large firing rate regime is 
captured by a series expansion of log() around 1. Standard asymptotic calculations enable equation (21) 
to be re-written as: 


Tmro{q,0) 


where 2 : 


- ]^{l + qxQ) 

(1 + qxoYe 

I2[qxi - d{l + ga;o)] 

Q 1 + 

qxi - 6»(1 + qxo) 


-0{z\l 


( 22 ) 


(23) 

(24) 


The modulation of the firing rate heterogeneity can be understood simply by the dominant term (i.e., 
first term): 

(25) 

specifically, the fraction q/0 is key because the second factor xi does not vary much as the correlation 
between q and 6 changes. 


^ Although ignoring the refractory period could be problematic for large firing rates, we emphasize that the value of our 
analysis is not in quantitative matching of simulations but rather for a deeper understanding of how network attributes 
effect the outputs. A similar calculation has been performed with the refractory period (not shown), but the asymptotic 
formulas are not insightful. 
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In the rhythmic firing regime, the modulation of the range of firing rates can be understood simply 
with range of values given by the fraction QjlOj. Specifically, the ratio qj^ yields values, and the 
range of these values for a particular parameter is indicative of the relative range of firing rates. In this 
regime (rhythmic firing) when ^(q, 0) < 0, the extreme values consists of: i) larger qj values that tend to 
occur with smaller Oj values, resulting in an amplification (and relatively larger) q^ for the upper range 
of values, and ii) smaller qj values and larger Oj^ resulting in an overall smaller (small times small) values 
to account for the lower range of q^ values. When p(q, 0) > 0, similar reasoning applies to larger qj 
values tending to occur with larger Oj values (large times small) and smaller values (small times large), 
resulting in a diminished range of values than when p > 0. Note that the mean of qj/Oj (across the 
Ne population) is approximately constant as g varies (not shown). 

A side-by-side comparison shows how similar the dynamics are (Fig. -B, E-F) and is validation 
that the range of excitatory firing rates is driven by this factor {q_/6). Figure [^compares the excitatory 
firing rate range of the LIF simulations (left column) to the analytic descriptions (right column, sections 


ae = 0.44 (dark orange). Of course we do not expect precise quantitative matching between the analytic 
description (Fig. right column) and the actual network simulation (Fig. left column) because the 
analysis had many assumptions meant to highlight a proof of principle. Nevertheless, the theory is 
quite descriptive and captures the general trend (the values of equation ( [^ are shown in Fig. [^,F 
for completeness). Therefore, in the rhythmic firing regime, q/0 compactly describes how when intrinsic 
and network heterogeneity are anti-correlated the range of firing rates is larger than when they are 
correlated. 


3.2 


3.3) as a function of p(q, 0) for two fixed levels of heterogeneity: aq = ae = I (black) and aq = 


3.3 Analytic description of heterogeneous firing rate range in asynchronous networks 


For the asynchronous regime, the range of firing rates actually increases as the correlation between 6 and 
q increases (Fig. [^), which is the opposite trend compared to when the network is firing rhythmically 
(Fig. 1^ and|^). This section provides an analytic description for this phenomena. 

In contrast to when the population firing rate is rhythmic, the asynchronous regime (i.e., power 
spectrum of population firing rate is fiat) has much less net synaptic input on average. Therefore, we 
cannot ignore the noise variable rfE and must consider a different regime than in the previous section. 
Here, the refractory period r^g/ is ignored, and all of the random variables are again assumed to be 
parameters except the noise variable rfE because of how crucial it is for firing in this regime. The formula 
for the reduced firing rate is (cf. equations 


/ .X l^qxo f 

W? 7 

where [ ] + represents thresholding: if rfE 
equation can be re-written by substituting 
integration: 


(16) and (^): 


n + 


log ( 


qxi-\-rrE-9(l-\-qxo) I 


qxi+rfE 


3 —E 


drfE 


(26) 


< 9{1 + qxo) — qxi it is 0, otherwise it is l/log(-). This 
[1/log(-)]+ with an integral, and interchanging the order of 


Tmro{q,9) 


1 +qxo 
(TeV^ 

1 + qxQ 

<xeV^ 


oo pqxi+rfE 


J I 

J —oo J q '. 

n 

J —oo J V 


'-oo J qxi-\-rrE-0(l-\-qxo) 
ny-qxi-\-9{l-\-qxo) 


-oo Jy-qxi 


M{y) dye drfE 

difEM{y) dy 


(27) 


where M{y) is an antiderivative of [l/log(-)] + . In the asynchronous firing regime, the modulation of 
the range of firing rates can be understood simply with the last term of upper limit of the inner integral: 


9{l + qxo) 


(28) 
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Fig. 5 Analytic description of excitatory firing rate ranges. Here, the firing rate ranges are for a fixed level of heterogeneity 
( 7 g = ag and different correlation ^(q, 0 ) values. A)-B) Noisy Rhythm regime: the range of firing rates for the completely 
heterogeneous (cTg = ctq = 1, black dots) network and with less heterogeneity (cTg = ag = 0.44, dark orange dots) are 
plotted as the correlation varies on the same axi s (se e left and right vertical axes for respective scales). B) The theory in 
section [3^ provides an analytic description (see ( |25| ) for how the correlation of intrinsic and network heterogeneity lead 
to relatively different firing rate ranges in this regime. As the correlation g increases, the range of firing rates tends to 
decrease. C)-D) Asynchrono us r egime: similar to A)-B), but the reduction theory in D) is in sectionagain providing 
an analytic description (see ( |28| >) for how the range of firing rates increases as the correlation g increases. E)-F) Sharp 
Rhythm regime: similar t o A) - B), using the same analytic description. The analytic descriptions in the right column (B, 
D, F) use xq and xi (see ( |19| >-( [20| ) values obtained from ^ = 0. Error bars i n A , C, and E are estimates of the standard 
deviation about the sample mean of the range (see last paragraph of Section l^for details). 
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because the other pieces of ro do not vary much across the neurons. Furthermore, since the focus is 
on understanding the range of firing rates as determined by the correlation of (q, 0), then the product: 
Oq is what determines the relative change in the range (like before, xq does not change much across 
the excitatory neurons). This is in sharp contrast to equation (25), where the intrinsic heterogeneity 
0 divides the network heterogeneity q. Figure [^-D shows that the analytic description (28) clearly 
captures the trend in the range of firing rates as g changes. 

These formulas ( [2^ , ( [2^ , clearly show how the two forms of heterogeneity (q, 9) effect one another 
to yield the different trends in the range of firing rates in different regimes. Simply put, in the rhythmic 
regime, the relationship q/0 determines how the range of the firing rates change while in the asynchronous 
regime, qO determines this; note that plotting q/0 and qO without xq and Xi does not appreciably change 
the shape of the curves in Figure(not shown). 

The key to our analysis was to focus on the main sources of heterogeneity, and therefore the firing 
rate heterogeneity, in a proper framework (probability density functions) to describe the essence of how 
the firing rate range changes. This analysis was successful partly because the heterogeneity (g, 0) drove 
the changes in the firing rate range, as opposed to other factors such as external noise, etc., and our 
analysis centered on these parameters. 


4 Discussion 


We studied how two forms of heterogeneity: intrinsic and network, effect the (excitatory) firing rate 
distribution of a recurrently coupled stochastic network of leaky integrate-and-fire neurons. Since the 
relationship between intrinsic and network heterogeneity is not known (to the best of our knowledge), 
we systematically varied the relationship or correlation to assess the effects on the network in different 
regimes. This work showed how the firing rate range changes with the correlation of intrinsic and 
network heterogeneity: in the rhythmic or oscillatory regime, the firing rate range tends to decrease with 
increasing correlation (i.e., when larger firing thresholds tend to have larger synaptic input amplification), 
while the opposite trend is observed in the asynchronous regime. These observations were captured by 


the analytic descriptions in equations (25) and (28). We also found that the firing rate ranges can be 


relatively large or small depending on the correlation between intrinsic and network heterogeneity, thus 
the overall level of heterogeneity can be mitigated or amplified depending on this relationship. If the 
relationship between intrinsic and network heterogeneity could be measured in a cortical network, these 
results would enable predictions for the range of response heterogeneity. Although we chose to analyze 
two specific forms of heterogeneity in a theoretical model, connections to experimental recordings of 
heterogeneity and firing rates may be possible even with related neural attributes (i.e., the membrane 
time constant as a proxy for firing threshold, or relating network connectivity to input variability). Also, 
the framework presented here could in principle be adapted to other heterogeneous neural attributes 
that would naturally require augmentations. 


Marder (2011) showed how combining intrinsic and synaptic conductance heterogeneity could lead 


to similar outputs, depending on their relationship. Their work was naturally different than the work 
here: different parameters were varied and the underlying neuron model had ionic currents and specific 
circuitry motivated by the crustacean stomatogastric ganglion. Also, they were interested in the rhythmic 
output of the network rather than characterizing the range of the population response heterogeneity. 


The results in Marder (2011) are similar in spirit to what has been shown here in that different sets of 
parameters can result in similar output (also see Marder and Goaillard ( 2QQ6| )); in our case, the range of 
excitatory firing rates can arise with different levels of heterogeneity by tuning the correlation between 
the two sources of heterogeneity. One of the conclusions of their work is that the intrinsic and network 
parameters, or heterogeneities, must be taken as a whole and the correlation among these parameters 
is crucial in determining network output. Our study compliments these assertions in a specific way, 
by determining how the correlation of two parameters alters the range of the excitatory firing rates in 
different regimes. 
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Unlike the work of Hermann and Touboul (2012); Mejias and Longtin (2012), we do not consider 


how synaptic heterogeneity can induce different dynamics but rather focus on specific spiking regimes 
that are similar with and without heterogeneities. Their work has a more detailed characterization of 
the dynamics, whereas our work explores the effects of two specific forms of heterogeneities away from 
the bifurcation points. Studying how the relationship of intrinsic and network heterogeneity induces 
different dynamics would be interesting but is beyond the scope of this paper. A recent paper of |Ostoji^ 
(2014) calls the first regime that we termed noisy rhythm ’a second type of asynchronous’ network but 
with stronger coupling. [Ostojic ( |2Q14 ) found that these two types of asynchronous networks (classic 
and strongly coupled) processed external stimuli differently, a feature that was not considered in this 
paper. These two regimes (Ostojic’s strongly asynchronous regime and the noisy rhythm here) are similar 
because the coupling is relatively strong and both autocorrelation functions of the population firing rates 
are similar. Consistent with illustrating that these regimes are different, we have shown how the classic 
asynchronous regime is different than the noisy rhythm or ’strongly coupled asynchronous’ (Ostojic) 
because the firing rate ranges change in distinct ways. 


A related study analyzes the interplay of two sources of intrinsic heterogeneity (Mejias and Longtin 
2Q14| ); specifically, Mejias and Longtin (2014) studied how heterogeneity in the excitatory and inhibitory 
spiking thresholds had different effects on a coupled network (LIT with excitatory and inhibitory neu¬ 
rons). They found distinct roles for heterogeneity in each type of neuron: excitatory heterogeneity can 
increase firing rate and linearizes output response, inhibitory heterogeneity can decrease network re¬ 


sponse and lead to gain control of input/output response (see Mejias and Longtin (2014) for further 
details). Their analysis also characterized the heterogeneity-induced transitions from asynchrony to syn¬ 
chrony and briefly considered the combined effects of these two attributes. Our work only considered 
excitatory heterogeneity; the effects of inhibitory neuron heterogeneity combined with network hetero¬ 
geneity are not known and a potential future direction of research. Another study by [Hnnsberger et al 


(2014) examined how varying both (white) noise in the voltage and heterogeneity in the threshold for 
firing led to different information (mutual information) content in spiking neuron models (LIT and 
Fitzhugh-Nagumo). They found an optimal level of heterogeneity for maximizing information content 
for a fixed level of noise (and for fixed level of heterogeneity there was an optimal level of noise, i.e., 
a stochastic resonance). Their results are distinct from [Mejias and Longtin (2012) and [Trip athy et al 


( 2Q13| ) (who also found optimal information by tuning parameters) because they considered the inter¬ 


play of those two sources of variability and determined that they interact nonlinearly (e.g., the optimal 
parameters are different with both components). In our study, we did not systematically analyze how 
varying the (correlated) noise level cfe and heterogeneity effect network statistics. [Hnnsberger et al] 
(2014) explained their simulation results by comparing how each component desynchronized and/or 
linearized the network response properties, whereas we provide an analytic explanation. Finally, [Leng^ 
et al ( 2013[ ) recently simulated an LIF network with a large number of heterogeneous intrinsic and net¬ 
work parameters. They find that heterogeneity can increase response time and paradoxically less variable 
responses (reliability), though they do not provide underlying mechanistic explanations for their results. 

Our study provides a more complete understanding of how heterogeneities interact and result in 
modulation of the firing rate statistics, which may ultimately lead to a better understanding of neural 
coding in coupled neural networks. Even though the firing rate is a first order measure of the response 
statistics, the range of this quantity has an impact on coding. There have been a number of recent studies 


focusing on the impact of heterogeneity on neural coding. Padmanabhan and Urban (2010) showed 


with recordings of mitral cells in mice olfactory bulb that heterogeneous cells had lower correlated 
activity, which is thought to increase information capacity of a given population. Similarly, [Chelarn and 


Dragoi (2008) found diverse orientation tuning curves enhances coding with increased information via 
a reduction in correlated activity of a coupled LIF network. Shamir and Sompolinsky (2006) proposed 


a theoretical explanation for the benefits of diversity/heterogeneity in population coding, whereby the 
information capacity is not limited to the correlation of activity. A future direction of study is how 
different forms of heterogeneities considered in this paper lead to changes in the second order statistics 
(correlation or co-variability), which also have implications for coding in neural systems. Although many 
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of the previous studies conclude that heterogeneity generally leads to lower co-variability (Chelaru and 


Dragoi 2008 Padmanabhan and Urban 2010), better discrimination (Marsat and Maler 


2010 Mejias 


et al 2013), and ultimately enhanced coding, the subtleties of how co-variability is modulated is not 


completely known and remains an active area of research (Ponce-Alvarez et al 2013 Ruff and Cohen 


2014 Mocholetal 2015) 
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Appendix: Generating correlated q and Q 

Given two vectors Q (intrinsic heterogeneity) and q (network heterogeneity), we can generate a new pair 
of vectors (of the same size) that have any desired correlation coefficient q G (—1,1). In this paper, we 
choose to keep q fixed and generate a new vector 'd that has the same sample mean (/i(^)) and sample 
standard deviation (cr(0)) of Q. Note that there are infinitely many ways to generate two such vectors 
if we only require that the mean and standard deviation of the new vectors be equal to the original 
statistics of the vector. The algorithm we use is as follows. 

— INPUTS: (q,6/,p) 

— Set Lp — cos“^(p) 

— Shift input vectors so they have zero mean: 

qo = q - /u(q) 

6>o = 6> - Ai(6>). 

— Calculate orthogonal complement to qo: 


— Create unit vectors out of qo and z: 
q = qo/l|qo||, z = z/||z|| 

— Create vector with prescribed correlation and zero mean: 
d = cos((/:?)q + sin((^)z 

— Sett9 = ^0 + MW 

— OUTPUT: where correlation coefficient of 'd and q is p, /i('i^) = /i(^), and (y(p&) = cr(0). 


References 

Apfaltrer F, Ly C, Tranchina D (2006) Population density methods for stochastic neurons with realistic 
synaptic kinetics: Firing rate dynamics and fast computational methods. Network: Computation in 
Neural Systems 17:373-418 

Borgers C, Kopell N (2003) Synchronization in networks of excitatory and inhibitory neurons with 
sparse, random connectivity. Neural computation 15(3):509-538 

Bremaud A, West D, Thomson A (2007) Binomial parameters differ across neocortical layers and with 
different classes of connections in adult rat and cat neocortex. Proceedings of the National Academy 
of Sciences 104:14,134-14,139 

Burton S, Ermentrout B, Urban N (2012) Intrinsic heterogeneity in oscillatory dynamics limits 
correlation-induced neural synchronization. Journal of Neurophysiology 108:2115-2133 

Buzsaki G, Wang XJ (2012) Mechanisms of gamma oscillations. Annual review of neuroscience 35:203 

Chelaru M, Dragoi V (2008) Efficient coding in heterogeneous neuronal populations. Proceedings of the 
National Academy of Sciences 105:16,344-16,349 

Chow CC (1998) Phase-locking in weakly heterogeneous neuronal networks. Physica D: Nonlinear Phe¬ 
nomena 118:343-370 


19 





































Economo MN, White JA (2012) Membrane properties and the balance between excitation and inhibition 
control gamma-frequency oscillations arising from feedback inhibition. PLoS computational biology 
8(l):el002,354 

Haskell E, Nykamp DQ, Tranchina D (2001) Population density methods for large-scale modeling of 
neuronal networks with realistic synaptic kinetics: cutting the dimension down to size. Network: 
Compututation in Neural Systems 12:141-174 

Hermann G, Touboul J (2012) Heterogeneous connections induce oscillations in large-scale networks. 
Physical Review Letters 109:018,702 

Hertag L, Durstewitz D, Brunei N (2014) Analytical approximations of the firing rate of an adaptive 
exponential integrate-and-fire neuron in the presence of synaptic noise. Erontiers in computational 
neuroscience 8 

Hunsberger E, Scott M, Eliasmith C (2014) The competing benefits of noise and heterogeneity in neural 
coding. Neural computation 26(8): 1600-1623 

Lengler J, Jug E, Steger A (2013) Reliable neuronal systems: The importance of heterogeneity. PloS one 
8(12):e80,694 

Levy RB, Reyes AD (2012) Spatial profile of excitatory and inhibitory synaptic connectivity in mouse 
primary auditory cortex. The Journal of Neuroscience 32(16):5609-5619 
Ly C (2013) A Principled Dimension Reduction Method for the Population Density Approach to Mod¬ 
eling Networks of Neurons with Synaptic Dynamics. Neural Computation 25:2682-2708 
Ly C (2014) Dynamics of coupled noisy neural oscillators with heterogeneous phase resetting curves. 

SIAM Journal on Applied Dynamical Systems 13:1733-1755 
Ly C, Tranchina D (2009) Spike Train Statistics and Dynamics with Synaptic Input from any Renewal 
Process: A Population Density Approach. Neural Computation 21:360-396, DOI 10.1162/neco.2008. 
03-08-743 

Ly C, Middleton J, Doiron B (2012) Cellular and circuit mechanisms maintain low spike co-variability 
and enhance population coding in somatosensory cortex. Erontiers in Computational Neuroscience 
6:1-26, DOI 10.3389/fncom.2012.00007 

Marder E (2011) Variability, compensation, and modulation in neurons and circuits. Proceedings of the 
National Academy of Sciences 108:15,542-15,548 

Marder E, Goaillard J (2006) Variability, compensation and homeostasis in neuron and network function. 
Nature Reviews Neuroscience 7:563-574 

Markram H, Liibke J, Erotscher M, Roth A, Sakmann B (1997) Physiology and anatomy of synaptic 
connections between thick tufted pyramidal neurones in the developing rat neocortex. The Journal of 
Physiology 500:409 

Marsat G, Maler L (2010) Neural heterogeneity and efficient population codes for communication signals. 
Journal of neurophysiology 104:2543-2555 

Mejias J, Longtin A (2012) Optimal heterogeneity for coding in spiking neural networks. Physical Review 
Letters 108:228,102 

Mejias J, Longtin A (2014) Differential effects of excitatory and inhibitory heterogeneity on the gain 
and asynchronous state of sparse cortical networks. Erontiers in computational neuroscience 8 
Mejias JE, Marsat G, Bol K, Maler L, Longtin A (2013) Learning contrast-invariant cancellation of 
redundant signals in neural systems. PLoS computational biology 9(9):el003,180 
Mochol G, Hermoso-Mendizabal A, Sakata S, Harris KD, de la Rocha J (2015) Stochastic transitions into 
silence cause noise correlations in cortical circuits. Proceedings of the National Academy of Sciences 
112(ll):3529-3534 

Moreno-Bote R, Parga N (2006) Auto- and Grosscorrelograms for the Spike Response of Leaky Integrate- 
and-Eire Neurons with Slow Synapses. Physical Review Letters 96:028,101 
Nesse WH, Borisyuk A, Bressloff P (2008) Eluctuation-driven rhythmogenesis in an excitatory neuronal 
network with slow adaptation. Journal of Gomputational Neuroscience 25:317-333 
Nicola W, Ly G, Gampbell SA (2015) One-dimensional population density approaches to recurrently 
coupled networks of neurons with noise. SIAM Journal on Applied Mathematics (in press) :- 


20 



Nykamp D, Tranchina D (2001) A Population Density Approach That Facilitates Large-Scale Modeling 
of Neural Networks: Extension to Slow Inhibitory Synapses. Neural Computation 13:511-546 

Ostojic S (2014) Two types of asynchronous activity in networks of excitatory and inhibitory spiking 
neurons. Nature neuroscience 17:594-600 

Ostojic S, Brunei N, Hakim V (2009) Synchronization properties of networks of electrically coupled neu¬ 
rons in the presence of noise and heterogeneities. Journal of Computational Neuroscience 26(3):369- 
392 

Oswald A, Doiron B, Rinzel J, Reyes A (2009) Spatial profile and differential recruitment of gabab 
modulate oscillatory activity in auditory cortex. The Journal of Neuroscience 29:10,321-10,334 

Padmanabhan K, Urban N (2010) Intrinsic biophysical diversity decorrelates neuronal firing while in¬ 
creasing information content. Nature Neuroscience 13:1276-1282 

Parker D (2003) Variable properties in a single class of excitatory spinal synapse. The Journal of neu¬ 
roscience 23(8):3154-3163 

Ponce-Alvarez A, Thiele A, Albright T, Stoner G, Deco G (2013) Stimulus-dependent variability and 
noise correlations in cortical mt neurons. Proceedings of the National Academy of Sciences 110:13,162- 
13,167 

Ruff DA, Gohen MR (2014) Attention can either increase or decrease spike count correlations in visual 
cortex. Nature neuroscience 17(11): 1591-1597 

Shamir M, Sompolinsky H (2006) Implications of neuronal diversity on population coding. Neural Gom- 
putation 18:1951-1986 

Strogatz SH, Mirollo RE (1991) Stability of incoherence in a population of coupled oscillators. Journal 
of Statistical Physics 63:613-635 

Tripathy S, Padmanabhan K, Gerkin R, Urban N (2013) Intermediate intrinsic diversity enhances neural 
population coding. Proceedings of the National Academy of Sciences 110:8248-8253 

Wang XJ (2010) Neurophysiological and computational principles of cortical rhythms in cognition. 
Physiological reviews 90:1195-1268 

Xue M, Atallah BV, Scanziani M (2014) Equalizing excitation-inhibition ratios across visual cortical 
neurons. Nature 511:596-600 

Yim M, Aertsen A, Rotter S (2013) Impact of intrinsic biophysical diversity on the activity of spiking 
neurons. Physical Review E 87:032,710 


21 



